{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-9-46044d6f247d>:49: UserWarning: In Matplotlib 3.3 individual lines on a stem plot will be added as a LineCollection instead of individual lines. This significantly improves the performance of a stem plot. To remove this warning and switch to the new behaviour, set the \"use_line_collection\" keyword argument to True.\n",
      "  ax[0].stem(freq_shift,amp_shifted)\n",
      "<ipython-input-9-46044d6f247d>:59: UserWarning: In Matplotlib 3.3 individual lines on a stem plot will be added as a LineCollection instead of individual lines. This significantly improves the performance of a stem plot. To remove this warning and switch to the new behaviour, set the \"use_line_collection\" keyword argument to True.\n",
      "  ax[1].stem(freq_shift,amp_shifted)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAAD4CAYAAACg9uHUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAeqElEQVR4nO3df5DkdX3n8ec7y5JakbvVMBAYSMBzmQqJp5tQEIpKakhyAlsmrF7iQVIRL9at5mJVsJQEohVNUndyofLjUhjIqhR4GoRKYCVx48oJc8QSRGAXFsQNCxrc2T22UAcZmdNled8f8x3tHbq3Z/v77fl+vz3PR1VXd39/9PfVM5/ve97f7v72RGYiSZKkevxQ3QEkSZJWMpsxSZKkGtmMSZIk1chmTJIkqUY2Y5IkSTU6qu4A3axduzZf/epX1x2jr+985zscc8wxdcfoy5zVa0vWtuR84IEHnsnMsbpzVMH6Va225IT2ZDVn9crWsEY2YyeccAL3339/3TH6mpqaYnJysu4YfZmzem3J2pacEfGvdWeoivWrWm3JCe3Jas7qla1hvk0pSZJUI5sxSZKkGtmMSZIk1chmTJIkqUY2Y5IkSTXq24xFxCkRcVdEPBYRj0bE7xbTXxkRd0TE48X1K3qsf0FE7IqI3RFxRdVP4Eht2T7NuVfdyWlXfJpzr7qTLdun644kNdIo7CujVr8kjaalvDL2AvDuzPwJ4GeB34mIM4ArgM9l5jrgc8X9Q0TEKuBDwIXAGcAlxbq12LJ9mitv3cn0zBwJTM/MceWtO1v5R0YaphHaV0amfkkaXX2bsczcl5kPFrefAx4DxoGLgBuLxW4ENnZZ/Sxgd2Y+mZnfAz5ZrFeLq7ftYu7AwUOmzR04yNXbdtWUSGqmUdlXRql+SRpdR/SlrxFxKrAe+CJwQmbug/mCFxHHd1llHPh6x/09wNk9HnsTsAlgbGyMqampI4m2JNMzcz2nD7K92dnZoeSsmjmr15asg+asel9pgrbXr6qN+hiuQ1uymrN5ltyMRcTLgb8HLsvMb0fEklbrMi27LZiZm4HNABMTEzmMb90dv/fOrn9kxteuGehbftvy7cDmrF5bsg6as+p9pW6jUL+qNupjuA5tyWrO5lnS2ZQRsZr5QvaJzLy1mPx0RJxYzD8R2N9l1T3AKR33Twb2Dh63nMvPn2DN6lWHTFuzehWXnz9RUyKpmUZpXxmV+iVpdC3lbMoAPgo8lpl/3jHrduDS4valwKe6rP4lYF1EnBYRRwMXF+vVYuP6cT74ptdw9Kr5pz2+dg0ffNNr2Lh+vK5IUiONyr4ySvVL0uhaytuU5wK/CeyMiB3FtD8ArgJuiYi3AU8BvwYQEScBH8nMDZn5QkS8E9gGrAKuz8xHq34SR2Lj+nFuuu8pAG5++zl1RpEabUT2lZGqX5JGU99mLDM/T/fPTgD8Ypfl9wIbOu5vBbYOGlCSBmX9ktQGfgO/JElSjWzGJEmSamQzJkmSVCObMUmSpBrZjEmSJNXIZkySJKlGNmOSJEk1shmTJEmqkc2YJElSjWzGJEmSamQzJkmSVCObMUmSpBrZjEmSJNXIZkySJKlGNmOSJEk1OqrfAhFxPfAGYH9m/lQx7WZgolhkLTCTma/rsu7XgOeAg8ALmXlmRbklSZJGQt9mDLgBuAb42MKEzPxPC7cj4s+AZw+z/nmZ+cygASWpDA8oJTVd32YsM++OiFO7zYuIAN4M/EK1sSSpMjfgAaWkBiv7mbGfA57OzMd7zE/gsxHxQERsKrktSTpimXk38M1u8zoOKG9a1lCS1GEpb1MeziUcvoidm5l7I+J44I6I+EpRGF+iaNY2AYyNjTE1NVUyWm8zM3MApbcxOzs71JxVMWf12pK1bM6q9pUGW+oBZQJ/k5mbuy20nPWrKitlDC+ntmQ1ZwNlZt8LcCrwyKJpRwFPAycv8TE+ALxnKcuefvrpOUxvvu4L+ebrvlD6ce66667yYZaBOavXlqxlc1a1r/QD3J9LqA2DXrrVsGL6tcC7D7PeScX18cBDwM/329aw61dVVsoYXk5tyWrO6pWtYWXepvwl4CuZuafbzIg4JiKOXbgNvB54pMT2JKkyEXEU8Cbg5l7LZObe4no/cBtw1vKkk7SS9G3GIuIm4B5gIiL2RMTbilkXs+gtyog4KSK2FndPAD4fEQ8B9wGfzszPVBddkkrxgFJSIyzlbMpLekx/a5dpe4ENxe0ngdeWzCdJpRQHlJPAcRGxB3h/Zn6UHgeUwEcycwPzB5S3zX/Gn6OAv/WAUtIwlP0AvyQ1mgeUkprOf4ckSZJUI5sxSZKkGtmMSZIk1chmTJIkqUY2Y5IkSTWyGZMkSaqRzZgkSVKNbMYkSZJqZDMmSZJUI5sxSZKkGtmMSZIk1chmTJIkqUY2Y5IkSTWyGZMkSaqRzZgkSVKN+jZjEXF9ROyPiEc6pn0gIqYjYkdx2dBj3QsiYldE7I6IK6oMLkmSNAqW8srYDcAFXab/RWa+rrhsXTwzIlYBHwIuBM4ALomIM8qElaQj5QGlpKbr24xl5t3ANwd47LOA3Zn5ZGZ+D/gkcNEAjyNJZdyAB5SSGuyoEuu+MyLeAtwPvDszv7Vo/jjw9Y77e4Czez1YRGwCNgGMjY0xNTVVItrhzczMAZTexuzs7FBzVsWc1WtL1rI5q9pX6pSZd0fEqQOs+v0DSoCIWDig/HJ16SRp8GbsWuBPgCyu/wz4rUXLRJf1stcDZuZmYDPAxMRETk5ODhitv2t33QPA5OQ5pR5namqKYeasijmr15asZXNWta80VGUHlMt5MFmVlXJAsZzaktWczTNQM5aZTy/cjogPA//YZbE9wCkd908G9g6yPUmqWKUHlMt5MFmVlXJAsZzaktWczTPQV1tExIkdd98IPNJlsS8B6yLitIg4GrgYuH2Q7UlSlTLz6cw8mJkvAh9m/i3JxTyglLQslvLVFjcB9wATEbEnIt4G/GlE7IyIh4HzgHcVy54UEVsBMvMF4J3ANuAx4JbMfHRIz0OSlswDSklN0vdtysy8pMvkj/ZYdi+woeP+VuAlZylJ0nIpDigngeMiYg/wfmAyIl7H/NuOXwPeXix7EvCRzNyQmS9ExMIB5Srgeg8oJQ1DmbMpJanxPKCU1HT+OyRJkqQa2YxJkiTVyGZMkiSpRjZjkiRJNbIZkyRJqpHNmCRJUo1sxiRJkmpkMyZJklQjmzFJkqQa2YxJkiTVyGZMkiSpRjZjkiRJNbIZkyRJqpHNmCRJUo1sxiRJkmrUtxmLiOsjYn9EPNIx7eqI+EpEPBwRt0XE2h7rfi0idkbEjoi4v8rgkrQU1jBJTbeUV8ZuAC5YNO0O4Kcy898D/wJceZj1z8vM12XmmYNFlKRSbsAaJqnB+jZjmXk38M1F0z6bmS8Ud+8FTh5CNkkqzRomqekiM/svFHEq8I+Z+VNd5v0DcHNmfrzLvK8C3wIS+JvM3HyYbWwCNgGMjY39zC233LLEp3DkPvjFOQCuPHtNqceZnZ3l5S9/eRWRhsqc1WtL1rI5q9pX+jnvvPMeGOYrT8OuYctZv6qyUsbwcmpLVnNWr3QNy8y+F+BU4JEu098L3EbR1HWZf1JxfTzwEPDzS9ne6aefnsP05uu+kG++7gulH+euu+4qH2YZmLN6bclaNmdV+0o/wP25hNow6GU5a9iw61dVVsoYXk5tyWrO6pWtYQOfTRkRlwJvAH6jCNKt0dtbXO8vCt5Zg25PkqpkDZPUFAM1YxFxAfD7wK9k5vM9ljkmIo5duA28Hnik27KStJysYZKaZClfbXETcA8wERF7IuJtwDXAscAdxSnf1xXLnhQRW4tVTwA+HxEPAfcBn87MzwzlWUhSD9YwSU13VL8FMvOSLpM/2mPZvcCG4vaTwGtLpZOkkqxhkprOb+CXJEmqkc2YJElSjWzGJEmSamQzJkmSVCObMUmSpBrZjEmSJNXIZkySJKlGNmOSJEk1shmTJEmqkc2YJElSjWzGJEmSamQzJkmSVCObMUmSpBrZjEmSJNXIZkySJKlGNmOSJEk16tuMRcT1EbE/Ih7pmPbKiLgjIh4vrl/RY90LImJXROyOiCuqDF6nLdunOfeqO3nrZ77DuVfdyZbt03VHkg7hGP0Ba9ihHBtqupU4RpfyytgNwAWLpl0BfC4z1wGfK+4fIiJWAR8CLgTOAC6JiDNKpW2ALdunufLWnUzPzAEwPTPHlbfuXBGDRe3gGH2JG7CGAY4NNd9KHaN9m7HMvBv45qLJFwE3FrdvBDZ2WfUsYHdmPpmZ3wM+WazXaldv28XcgYOHTJs7cJCrt+2qKZF0KMfooaxhP+DYUNOt1DF61IDrnZCZ+wAyc19EHN9lmXHg6x339wBn93rAiNgEbAIYGxtjampqwGj9zRQd9yDbWOjWu00fZuYyZmdnG5utU1tyQrOzVjlGy+wrDVdpDVvO+lWG9Wu42pK1yTnbOEarMGgzthTRZVr2WjgzNwObASYmJnJycnJIseDaXfcAMDl5zhGvO37vnV0Hy/jaNQwzcxlTU1ONzdapLTmh2VmrHKNl9pURsOQatpz1qwzr13C1JWuTc7ZxjFZh0LMpn46IEwGK6/1dltkDnNJx/2Rg74Dba4zLz59gzepVh0xbs3oVl58/UVMi6VCO0SVZkTXMsaGmW6ljdNBm7Hbg0uL2pcCnuizzJWBdRJwWEUcDFxfrtdrG9eN88E2v4ehV8z+68bVr+OCbXsPG9eM1J5PmOUaXZEXWMMeGmm6ljtG+b1NGxE3AJHBcROwB3g9cBdwSEW8DngJ+rVj2JOAjmbkhM1+IiHcC24BVwPWZ+ehwnsby2rh+nJvue4qZmRm2/f4v1B1HegnH6A9Yww7l2FDTrcQx2rcZy8xLesz6xS7L7gU2dNzfCmwdOJ0klWQNk9R0fgO/JElSjWzGJEmSamQzJkmSVCObMUmSpBrZjEmSJNXIZkySJKlGNmOSJEk1shmTJEmqkc2YJElSjWzGJEmSamQzJkmSVCObMUmSpBrZjEmSJNXIZkySJKlGNmOSJEk1GrgZi4iJiNjRcfl2RFy2aJnJiHi2Y5k/LB9ZkiRpdBw16IqZuQt4HUBErAKmgdu6LPrPmfmGQbcjScMQERPAzR2TXgX8YWb+Zccyk8CngK8Wk27NzD9etpCSVoSBm7FFfhF4IjP/taLHk6Sh8oBSUlNU9Zmxi4Gbesw7JyIeioh/ioifrGh7klQlDygl1ab0K2MRcTTwK8CVXWY/CPx4Zs5GxAZgC7Cux+NsAjYBjI2NMTU1VTZaTzMzcwCltjEzM8fBgweHmrMqs7Oz5qxYG7JWMUar2Fdaou8BJbAXeE9mPrp4geWsX1Wwfg1HW7K2IWebxmgVqnib8kLgwcx8evGMzPx2x+2tEfHXEXFcZj7TZdnNwGaAiYmJnJycrCBad9fuugeAyclzSj3GzMwMw8xZlampKXNWrA1ZqxijVewrTVfFAeVy1q8qWL+Goy1Z25CzTWO0ClW8TXkJPY4oI+JHIyKK22cV2/tGBduUpKoc9oAyM2eL21uB1RFx3HIHlDTaSr0yFhEvA/4D8PaOae8AyMzrgF8FfjsiXgDmgIszM8tsU5IqdtgDSuDpzEwPKCUNS6lmLDOfB35k0bTrOm5fA1xTZhuSNCweUEpqgqq+2kKSWscDSklN4L9DkiRJqpHNmCRJUo1sxiRJkmpkMyZJklQjmzFJkqQa2YxJkiTVyGZMkiSpRjZjkiRJNbIZkyRJqpHNmCRJUo1sxiRJkmpkMyZJklQjmzFJkqQa2YxJkiTVyGZMkiSpRqWasYj4WkTsjIgdEXF/l/kREX8VEbsj4uGI+Oky25OkKlnDJDXBURU8xnmZ+UyPeRcC64rL2cC1xbUkNYU1TFKthv025UXAx3LevcDaiDhxyNuUpKpYwyQNXdlXxhL4bEQk8DeZuXnR/HHg6x339xTT9i1+oIjYBGwCGBsbY2pqqmS03mZm5gBKbWNmZo6DBw8ONWdVZmdnzVmxNmStYoxWsa80XCU1bDnrVxWsX8PRlqxtyNmmMVqFss3YuZm5NyKOB+6IiK9k5t0d86PLOtntgYoiuBlgYmIiJycnS0br7dpd9wAwOXlOqceYmZlhmDmrMjU1Zc6KtSFrFWO0in2l4SqpYctZv6pg/RqOtmRtQ842jdEqlHqbMjP3Ftf7gduAsxYtsgc4peP+ycDeMtuUpKpYwyQ1wcDNWEQcExHHLtwGXg88smix24G3FGck/SzwbGa+5C1KSVpu1jBJTVHmbcoTgNsiYuFx/jYzPxMR7wDIzOuArcAGYDfwPPCfy8WVpMpYwyQ1wsDNWGY+Cby2y/TrOm4n8DuDbkOShsUaJqkp/AZ+SZKkGtmMSZIk1chmTJIkqUY2Y5IkSTWyGZMkSaqRzZgkSVKNbMYkSZJqZDMmSZJUI5sxSZKkGtmMSZIk1chmTJIkqUZl/lG4NJAt26e5etsu9s7McdLaNVx+/gQb148f8TKSVId+9cn6pSNlM6ZltWX7NFfeupO5AwcBmJ6Z48pbdwKwdgnLWNAk1alfDbN+aRC+TalldfW2Xd8vUgvmDhzk6m27jmgZSapDv/pk/dIgbMa0rPbOzPWdvpRlJKkO/eqT9UuDsBnTsjpp7Zq+05eyjCTVoV99sn5pEAM3YxFxSkTcFRGPRcSjEfG7XZaZjIhnI2JHcfnDcnHVdpefP8Ga1asOmbZm9SouP3/iiJaRyrKGaRD96pP1S4Mo8wH+F4B3Z+aDEXEs8EBE3JGZX1603D9n5htKbEcjZOEDrL/3dw/zvYMvMt5xptHU1ON9l5EqZA3TEetXw6xfGsTAzVhm7gP2Fbefi4jHgHFgcSGTDrFx/Tg33fcUADe//ZyBl5HKsIZpUP3qk/VLR6qSr7aIiFOB9cAXu8w+JyIeAvYC78nMR3s8xiZgE8DY2BhTU1NVROtqpvggZZltzMzMcfDgwaHmrMrs7Gzjcnb7HSzOWcXvaVia+DNdrIox2uTfQZXK1rDlrF9VsH6V16+GNXnfaerPtFObxmgVSjdjEfFy4O+ByzLz24tmPwj8eGbORsQGYAuwrtvjZOZmYDPAxMRETk5Olo3W07W77gFgcnLwI5Zrd93DzMwMw8xZlampqcbl7PY7WJyzit/TsDTxZ7pYFWO0yb+DqlRRw5azflXB+lVevxrW5H2nqT/TTm0ao1UodTZlRKxmvoh9IjNvXTw/M7+dmbPF7a3A6og4rsw2Jakq1jBJTVDmbMoAPgo8lpl/3mOZHy2WIyLOKrb3jUG3KUlVsYZJaooyb1OeC/wmsDMidhTT/gD4MYDMvA74VeC3I+IFYA64ODOzxDYlqSrWMEmNUOZsys8D0WeZa4BrBt2GJA2LNUxSU/gN/JIkSTWq5KstNLq2bJ/m6m272Dszx0kt//LCUXoukvobpX1+lJ6LXspmTD1t2T7NlbfuZO7AQQCmZ+a48tadAK0rAqP0XCT1N0r7/Cg9F3Xn25Tq6eptu76/8y+YO3CQq7ftqinR4EbpuUjqb5T2+VF6LurOZkw97S2+QXqp05tslJ6LpP5GaZ8fpeei7mzG1NNJa9cc0fQmG6XnIqm/UdrnR+m5qDubMfV0+fkTrFm96pBpa1av4vLzJ2pKNLhRei6S+hulfX6Unou68wP86mnhg6G/93cP872DLzLe4jN4Rum5SOpvlPb5UXou6m5FNWMLpwZPF++zn3rFpx3UfWxcP85N9z0FwM1vb94/vD0So/Rchs19RaNglPb5UXoueqkV04wtPjV4gacIS4dyX5Gk5bViPjPW7dTgBZ4iLP2A+4okLa8V04z1OwXYU4Slee4rkrS8VkwztvZlqw87/9+uOfz8BVu2T3PuVXfyxa9+kyeefZEt26eriCdV7pnnvssTz77IaVd8mnOvunPJY7WqfUWStDQrohl735adfOv5A4dd5tv/70DfP1YLn6VZ+FDzCy/CZTfvYP0ff9amTI2xZfs0r/ujz/LEM9/hhRch+cHnvfqN06r2FUnS0o18M/a+LTv5+L1P9V3uxYR33bzjsH9k/ugfHu36WZpvPX+Ad928g/dt2Vkqq1TW+7bs5F0372Bm7qUN1dyBg3zg9kd7rvsbH76nsn1FkrR0rTybcv4VqoeZO/BipY+bzL/SddnNOwZa9+P3PvWSP2bHHL2K//bG13j2mSqzZfs0H7j90a4NVz8zcwc49YpPl85wJPvKK162mvf/8k+6D0hSD5GZg68ccQHwP4FVwEcy86pF86OYvwF4HnhrZj7Y73F/+MR1eeKlfzlwLkntsu/Gy/juvsdjubc7jBpm/ZJWnrI1bOC3KSNiFfAh4ELgDOCSiDhj0WIXAuuKyybg2kG3J0lVsoZJaooyb1OeBezOzCcBIuKTwEXAlzuWuQj4WM6//HZvRKyNiBMzc9/hHvjk5/bz3//5r0tEk9Qmb6lns0OpYdYvaeUpW8PKfIB/HPh6x/09xbQjXQaAiNgUEfdHxP0lMi1ZLPsbIlL7RMCq0d1XKqthy12/JI2WMq+MdSvRiz+AtpRl5idmbgY2w/xnLn7/5/5riWi9dX6gflgnAkij4Nx/90o+8V/m/wfe0PeVJy8bzuMeXmU1bLnql6SGKlnDyjRje4BTOu6fDOwdYJll0e2Mro3rx7/flB3J2WkLDd2XH/sytzyeA53VJg3TDwX8+tk/xsue/79HPEar3FcarlU1TNLoKtOMfQlYFxGnAdPAxcCvL1rmduCdxWcxzgae7fd5sSNR1SnzC39ojtTaZx/nD359suu8EfujpYbqtw9MTX2j5xgdxJHsKy3YB2qvYZIEQGYOfGH+dO9/AZ4A3ltMewfwjuJ2MH+20hPATuDMpTzu6aefnm1w11131R1hScxZvbZkbUtO4P4sUYsGvQyjhlm/qtWWnJntyWrO6pWtYaW+9DUztwJbF027ruN2Ar9TZhuSNCzWMElNMPL/DkmSJKnJbMYkSZJqZDMmSZJUI5sxSZKkGpX6R+HDEhHPAbvqzrEExwHP1B1iCcxZvbZkbUvOicw8tu4QVbB+Va4tOaE9Wc1ZvVI1rNTZlEO0KzPPrDtEPxFxvzmr05ac0J6sbcpZd4YKWb8q1Jac0J6s5qxe2Rrm25SSJEk1shmTJEmqUVObsc11B1gic1arLTmhPVnNufza8lzMWb22ZDVn9UplbeQH+CVJklaKpr4yJkmStCLYjEmSJNWo1mYsIn4tIh6NiBcj4syO6adGxFxE7Cgu13XM+5mI2BkRuyPiryIi6spZzLuyyLIrIs6vM2eX3B+IiOmOn+OGfrnrEhEXFFl2R8QVdefpFBFfK36XOxZOX46IV0bEHRHxeHH9ipqyXR8R+yPikY5pPbPV9XvvkbM147ObttSvw2Ut5jWyhrVpfFi/Bs5m/VqQmbVdgJ8AJoAp4MyO6acCj/RY5z7gHCCAfwIurDHnGcBDwA8DpwFPAKvqytkl9weA93SZ3jN3TeNgVZHhVcDRRbYz6hybi/J9DThu0bQ/Ba4obl8B/I+asv088NOd+0uvbHX+3nvkbMX4PMxzakX96pO1sTWsLePD+lUqm/WruNT6ylhmPpaZS/6m6og4Efg3mXlPzj/rjwEbhxawcJicFwGfzMzvZuZXgd3AWXXlPAJdc9eY5yxgd2Y+mZnfAz5ZZGyyi4Abi9s3UtPvNzPvBr65aHKvbLX93nvk7KVp47OrttQvGLka1rTxYf0akPXrB5r8mbHTImJ7RPyfiPi5Yto4sKdjmT3FtLqMA1/vuL+Qp0k53xkRDxcvsy683Nsrd12almexBD4bEQ9ExKZi2gmZuQ+guD6+tnQv1StbE3/ObRifg2hD/YLm17A2jI+m5VnM+jU8lY3Pof87pIj438CPdpn13sz8VI/V9gE/lpnfiIifAbZExE8y/3L5YpV8N8eAOXvlGVrOlwQ4TG7gWuBPim3/CfBnwG8tZ74lalqexc7NzL0RcTxwR0R8pe5AA2raz7nx47Mt9QvaWcOsX8vC+jUclY7PoTdjmflLA6zzXeC7xe0HIuIJ4HTmO8yTOxY9GdhbV84izyld8gwt52JLzR0RHwb+sbjbK3ddmpbnEJm5t7jeHxG3Mf+S89MRcWJm7ive0tlfa8hD9crWqJ9zZj69cLup47Mt9avYVutqmPVr+Kxfw1F1/Wrk25QRMRYRq4rbrwLWAU8WL1k+FxE/W5zZ8xag1xHfcrgduDgifjgiTity3teUnMVAXvBGYOFMkK65lztfhy8B6yLitIg4Gri4yFi7iDgmIo5duA28nvmf4+3ApcVil1LvOFysV7ZG/d5bND6PSIvqFzS4hrVofFi/qrUy69dynInQ61I8gT3MH0U+DWwrpv9H4FHmz0h4EPjljnXOLJ70E8A1FP9FoI6cxbz3Fll20XG2UR05u+T+X8BO4OFigJzYL3eNY2ED8C9FpvfWnacj16uKcfhQMSbfW0z/EeBzwOPF9StryncT82+LHSjG6NsOl62u33uPnK0Znz2eUyvq1+GyHu5nXXcNa9P4sH4NnM/6VVz8d0iSJEk1auTblJIkSSuFzZgkSVKNbMYkSZJqZDMmSZJUI5sxSZKkGtmMSZIk1chmTJIkqUb/H94gBd+YcPBbAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np  # 导入 numpy 库，用于进行数学计算\n",
    "import matplotlib.pyplot as plt  # 导入 matplotlib.pyplot 库，用于绘制图形\n",
    "import math  # 导入 math 库，用于数学计算\n",
    "\n",
    "# 定义一个函数，用于生成正弦波\n",
    "def getSin(amp, freq, phase, sampleList):\n",
    "    return amp * np.sin(2 * math.pi * freq * sampleList + phase)\n",
    "# 定义一个函数，用于生成余弦波\n",
    "def getCos(amp, freq, phase, sampleList):\n",
    "    return amp * np.cos(2 * math.pi * freq * sampleList + phase)\n",
    "def denoise(arr, thresh):\n",
    "    mask= (arr > thresh)\n",
    "    mask= mask + 0 # 用mask加0把True和False转化成1和0显示\n",
    "    # 掩模的值是1，相乘保留数值。反之是0，相乘则为0\n",
    "    return mask * arr\n",
    "srate = 3000  # 设置采样率为 1000\n",
    "t = np.linspace(0, 1, srate)  # 在 0 到 1 的范围内生成等间距的采样点\n",
    "\n",
    "# 分别生成四个正余弦波\n",
    "s1 = getSin(amp=1.5, freq=30, phase=0, sampleList=t)\n",
    "s2 = getCos(amp=3, freq=5, phase=0, sampleList=t)\n",
    "s3 = getSin(amp=10, freq=100, phase=0, sampleList=t)\n",
    "s4 = getCos(amp=20, freq=120, phase=0, sampleList=t)\n",
    "m = s1 + s2 + s3 +s4 # 将四个正余弦波相加，得到混合波形\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# 获得傅里叶系数\n",
    "# 第一个参数是混合波形，第二个参数是“飞轮”测试频率上限（即采样率）\n",
    "fCoefs = np.fft.fft(m,srate)\n",
    "\n",
    "# 获得振幅列表：每一个绕线的重心到原点的距离\n",
    "amp_list=2 * np.abs(fCoefs/srate)\n",
    "\n",
    "# 把频率轴变换\n",
    "freqs = np.fft.fftfreq(len(amp_list), 1/srate)\n",
    "\n",
    "# 然后把 频率轴 和 数据 都变成 0hz 在中间，向左是负频率，向右是正频率的形式\n",
    "amp_shifted=np.fft.fftshift(amp_list)\n",
    "freq_shift=np.fft.fftshift(freqs)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "#去噪前\n",
    "fg ,ax = plt.subplots(1,2,figsize=(10,4))\n",
    "\n",
    "ax[0].stem(freq_shift,amp_shifted)\n",
    "ax[0].set_xlim([-150,150])\n",
    "ax[0].grid()  # 添加网格线\n",
    "\n",
    "\n",
    "#去噪\n",
    "amp_shifted[(freq_shift>110) | (freq_shift<-110)]= 0  #去除频率大于110hz和小于-110hz的噪音\n",
    "amp_shifted = denoise(amp_shifted, 1) #将振幅小于1的噪声全部去除   \n",
    "\n",
    "#去噪后\n",
    "ax[1].stem(freq_shift,amp_shifted)\n",
    "ax[1].set_xlim([-150,150])\n",
    "ax[1].set_ylim([-1,21])\n",
    "ax[1].grid()  # 添加网格线\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
